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Abstract. A new approach is proposed for reconstruction of images from 
Radon projections. Based on Fourier expansions in orthogonal polynomials 
of two and three variables, instead of Fourier transforms, the approach pro- 
vides a new algorithm for the computed tomography. The convergence of the 
algorithm is established under mild assumptions. 



The fundamental problem of the computed tomography (CT) is to reconstruct 
an image from its Radon projections (X-rays). To be more precise, let us consider 
the case of a two dimensional image, described by a density function f{x, y) defined 
on the unit disk — {(a;, y) : + y^ < 1} of the plane. A Radon projection 
of / is a line integral, 



where l{0,t) — {{x,y) : x cos + y sin 6 — t} O is a. line segment inside B^ . The 
reconstruction problem of CT requires finding / from its Radon projections. This 
question is completely solved if the complete data is given, that is, if projections for 
all t and 9 are given. In practice, however, only a finite number of projections can be 
measured. Hence, the essential problem of CT is to find an effective algorithm that 
produces a good approximation to / based on a finite number of Radon projections. 
The arrangement of the projections is referred as scanning geometry since it is 
determined by the design of the scanner. 

Currently, the most important method in CT is the filtered backprojection (FBP) 
method which, like several other methods, is based on techniques of Fourier trans- 
forms. In principle, the FBP method works for banded limited functions. It requires 
choosing a function Wb (a low-pass filter with cut-off frequency b) that approximates 
the (5-distribution, and it includes steps of linear interpolation and a discrete con- 
volution that approximates a continuous convolution. Each of these steps could 
introduce serious errors in the algorithm. See the discussion in 01^1 ■ 

The purpose of the present paper is to develop a direct approach for reconstruc- 
tion of images in CT. Instead of using Fourier transforms, we will use orthogonal 
expansions based on orthogonal polynomials on B^. Let 5„/ denote the n-th partial 
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1. Introduction 
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sum of such an expansion. It turns out that S2mf satisfies the following remarkable 
formula 

2m ^ „i 

where (j)^, — 2i/7r/(2m + 1) are equally spaced angles along the circumference of 
the disk and are polynomials of two variables given by explicit formulas (see 
Theorem 13.3(1 . Applying appropriate quadrature formulas to the integrals in ((l.l|l 
leads to an approximation, A2mf, of / that uses discrete Radon data. For example, 
using an appropriate Gaussian quadrature formula leads to an A2mf that uses the 
Radon data 

(1.2) {TZ^Jf-cos^): Q<v<2m, l<]<2m}, 

of the scanning geometry of parallel beams, and the polynomial A2mf is given by 
an explicit formula 

2m 2m 

(1.3) A2m{f; x,y)^Y.T. ^^^(/^ 2^)TiAx, y). 

u=Oj = l 

where Tj_^{x,y) are polynomials that are given by simple formulas. The operator 
A2mf provides a direct approximation to /, which is the essence of our reconstruc- 
tion algorithm. The new algorithm can be implemented easily as the Tj ^, are fixed 
polynomials that can be stored in a table beforehand. Furthermore, the construc- 
tion allows us to add a multiplier factor and there is a natural extension of the 
algorithm to a cylindrical domain in R"^. 

The operator .42m in (|1.3|) reproduces polynomials of degree 2m — 1. Fur- 
thermore, we will prove that the operator norm of ■42m in the uniform norm is 
0{m\og{m + 1)), only slightly worse than the norm of S'2m, which is 0{m). As a 
consequence, it follows that the algorithm converges uniformly on if / is a 
function. There is no need to assume that / is banded limited. 

At the moment the author is working with Dr. Christoph Hoeschen and Dr. Oleg 
Tischenko, Medical Physics Group of the National Research Center for Environment 
and Health, Germany, to implement the algorithm numerically. The results that 
we have obtained so far are very promising and will be reported elsewhere. 

The paper is organized as follows. The background of Radon transforms and 
orthogonal polynomials is discussed in the next section, where an identity that will 
play a fundamental role in the analysis is also proved. The Fourier orthogonal 
expansions are studied in the Section 3, where the identity (|1.1|) will be proved. 
The reconstruction algorithms are presented in Section 4. The convergence of the 
basic algorithm for the 2-D reconstruction is proved in Section 5 . 

2. Radon transform and polynomials 

Let = {{x,y) : x^ + y^ < 1} denote the unit disk on the plane. It is often 
convenient to use the polar coordinates 

a: = rcos6', y = rsm9, r>0, 0<6<27t. 

Let 9 be an angle as in the polar coordinates; that is, 9 is measured counterclockwise 
from the positive x-axis. Let i denote the line £{9,t) — {{x,y) : xcos9 + ysm9 = t} 
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for — 1 < i < 1. Clearly the line is perpendicular to the direction (cos 0, sin 6*) and 
\t\ is the distance between the line and the origin. We will use the notation 

(2.1) i{e,t) = e{e,t)nB^, o<e<2'K, -i<t<i, 

to denote the line segment of i inside B"^. The points on I{6,t) can be represented 
as follows: 

X = t cos ^ — s sin ^, y = tsinO + s cos 9, 

for s e [~VT^, vT^^ . 

The Radon projection of a function / in the direction 6 with parameter i e [—1,1] 
is denoted by TZ0{f; t), 

(2.2) neif;t):^f f{x,y)dxdy 

Ji(e,t) 



J ^ f{tcos6 — ssm6,tsm0 + scos6)ds. 



In the literature a Radon projection of a bivariate function is also called an X-ray. 
The definition shows that TZeif; t) = 7?.^+6i(/; —t). 

Let denote the space of polynomials of two variables and let 11^^ denote the 
subspace of polynomials of total degree n in 11^ , which has dimension dim 11^ = 
(n + l)(n + 2)/2. If P e then 

n k 
fc=0 j=0 

Let p be a polynomial in one variable and 9 G [0,7r]. For ^ = (cos 6*, sin 6*) and 
X = {x,y), the ridge polynomial p{6; •) is defined by 

p{6; x, y) := p((x, ^)) = p{x cos9 + y sin 9). 

This is a polynomial of two variables and is in 11^ if p is of degree n. It is constant 

on lines that are perpendicular to the direction (cos 0, sin 0). Especially important 
to us are the polynomials Uk{9;x,y), where Uk denotes the Chebyshev polynomial 
of the second kind, 

sin(fc + 1)9 

Uk{x) = ^-^ , X = cos9. 

sm9 

The polynomials Uk, k = 0,1,2, .. ., are orthogonal polynomials with respect to the 
weight function y/l — x-^ on [—1, 1]: 

— / Uk{x)Uj{x)\^ 1 — x'^dx = Sk j, k,j>0. 

J-1 

Ridge polynomials play an important role in our study of the Radon transform, as 
can be seen from the following simple fact. 

Lemma 2.1. For f G L^iB^), 

(2.3) / f{x,y)Uk{^;x,y)dxdy= f n^{f;t)Uk{t)dt. 

Jb^ J-1 
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Proof. The change of variables t — xcos0 + ysiii<p and s = — xsmt/f + ycosc 
amounts to a rotation, which leads to 



f{x,y)Uk{4>]x,y)dxdy = I f {t cos (j) — s siri 4>,t syo. (j) + s cos 4>)Uk{t)dtds 

f(t cos (j) ~ s sin0, f sin 4> + s cos (f>)dsUk(t)dt, 

the inner integral is exactly TZ^{f\ t) by (|2.2|l . □ 

Let Vk{B^) denote the space of orthogonal polynomials of degree k on with 
respect to the unit weight function; that is, P G Vk{B^) if P is of degree k and 

/ P{x, y)Q{x, y)dxdy = 0, for all Q e l^i^i- 
Jb^ 

A set of polynomials {Pj.k : < j < fc} in Vk{B^) is an orthonormal basis oiVk{B^) 
if 

-/ Pi.k{x,y)Pj,k{x,y)dxdy ^ 5i,j, 0<i,j<k. 
Jb^ 

It is known that the polynomials in Vk{B^) are eigenfunctions of a second order 
differential operator T) (see, for example, 

VP = -{k + 2){k-l)P, for all P eViB"^)] 

the operator V is defined by 

(2.4) A2-(x,V)2-2(x,V), 

where x — {x,y), V = {di,d2) is the gradient operator with di = d/dx and 
82 = d/dy, and A = 9^ + 9| is the usual Laplacian. It turns out that the Radon 
projections of orthogonal polynomials in V„(i3^) can be computed explicitly. 

Lemma 2.2. If P E Vk{B'^) then for each t e (-1, 1), < 6* < 2tt, 

2 

(2.5) ne(P]t) = y/l - t^Uk{t)P{cose,sme). 

K + 1 

Proof. A change of variable in H2.2|l shows that 
TZe{P;t) = \J\ - t2 y" P [tcosd- s\J\ ~ sin 6*, i sin 6* + sy/l - cos 9^ ds. 



The integral is a polynomial in t since an odd power of y/l — t in the integrant 
is always companied by an odd power of s, which has integral zero. Therefore, 
Q{t) := TZe{P; t) / \/l — t'^ is a polynomial of degree k int for every 0. Furthermore, 
the integral also shows that Q(l) = P(cos6', sin6'). By H2.3() . 

1 TMP^-^ 

-1 VT^ 



Uj (t) \J\ - t^dt = / P{x, y)Uj {9; x, y)dxdy = 0, 



52 



for j = 0, 1, . . . , A: — 1, since P G Vfc(i?^). Since Q is of degree /c, we conclude that 
Q{t) — cUk{t) for some constant independent of t. Setting t = \ and using the fact 
that Uk{l) = fc + 1, we have c = P(cos6',sin6')/(fc + 1). □ 
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This lemma was proved in It plays an important role in our development 
below. We have included its short proof in ordc to make the paper self-contained. 

There are several orthogonal or orthonormal bases that are known explicitly for 
Vfe(-B^) (see 15]). We will work with a special orthonormal basis that is given in 
terms of the ridge polynomials defined above. Setting f{x,y) = Uk{0]x,y) in (|2.3() 
and using H2.5(l . we derive from the orthogonality of the Chebyshev polynomials 
that 

(2.6) -( Uk{0;x,y)Uk{<P;x,y)dxdy^-^Uk{cos{<f>^9)). 

TT Jb2 k+l 

Recall that the zeros of Uk are cosOj^k, 1 < j < A:, where 9j^k = JT^/ik + 1). As a 
consequence of H2.6f) . we have the following result see also 

Proposition 2.3. An orthonormal basis ofVk{B^) is given by 



h ■■= {Uk{e.jy.x,y) : < J < fc}, 0j,fe 



fc + 1' 

In particular, the set {P^ : < k < n} is an orthonormal basis for 11^^. 

The polynomials Uk{4>]x,y) also satisfy a discrete orthogonality. Let 

2tiv 

Q<v <2m. 



2m + 1 

The discrete orthogonality will follow from the following identity: 
Proposition 2.4. For fc > and 9 G [0, 27r], 

2m 



(2.7) —^—y^Uki(l)^;cos9,sm9)Uk{(l>^;x,y)^Uki9;x,y), x,yeB^. 
Proof. In order to prove (|2.7|l we will need the elementary identities 

n n 

(2.8) ^^sin/c0j, — and ^^cos/c^^ = 



0, if fc 7^ mod n + 1 

n + 1 , if fc = mod n + 1 



that hold for all nonnegative integers fc. 

Let us denote by I}^ the left hand side of H2.7II . First we consider the case of 
k = 21. Since U21 is an even polynomial, we can write it as U2i {t) — X]j"=o ^j^^"' - 
Using the polar coordinates x = r cos cj) and y — rsin0, and the fact that (cos^')^"' 
can be written as a sum of cos2i?/', we can rearrange the sum to get 

/ 

(2.9) U2i{9;x,y) = U2i{rcos{9 - 0)) = ^ 6,r2J(cos(0 - cj))f^ 

3=0 
I 

= Y,h,{r)cos2j{9-4,), 

where bj(r) is a polynomial of degree 2j in r. Furthermore, we have 

/ 

(2.10) C/2;('/'^;cos6',sin6») = U2i{cos{(j)^ ~ 9)) = 1 + 2^cos2j{9 - (j)^). 
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These two expressions allow us to write 

II ^ 2m 

hi = ^di^6j(r)^—-— ^ cos 21(6* cos 2j(0- 

i=0 j"=0 i^=t) 

where do = 1 and di = 2 for I < i < I. The addition formula of the cosine function 
shows that 

2 cos 2i{9 — (f)^) cos 2j(0 — 01,) 

= cos2((i + j)0„ - iO ~ + cos2((i - j)(/)„ - i9 + j4>) 
= cos2(i + cos2{i9 + jcf)) + sin2(i + j)(f)i, s\n2[i9 + j(f)) 
+ cos2(i — cos2{i9 — jcf)) + sin2(i — sin 2[i9 — jcf)). 

Since n + 1 = 2m + 1 is odd and cos2(i ± i)4>u have even indices, it follows from 
the above elementary trigonometric identity and (|2.8() that 



^ 2m I 0, if«7^j, 

2j„TT^^°'^*^^~'^"^''°''^^^'^~'^"^ " 1 5Cos2j(0-0), ifi^j^O, 



"^=0 If, ifi=j = 0. 

Consequently, using H2.9|l again, we conclude that 
I 

hi = W co*^ -</>) = t/2i(^ cos(0 - 0)) = c/2i(^; a;, y)- 

This completes the proof for the case k — 21. For the case fc = 2Z — 1, we need to 
use the identities 

I 

U2i-i{(f>u;x,y) = U2i{rcos{(f>^ - (/))) = ^5j(r)cos(2j - - cj)), 

where bj{r) is a polynomial of degree 2j — 1 in r, which is derived using the fact 
that (cos6')^''~^ can be written as a sum of cos(2j — 1)9. Furthermore, we have 

i72/_i cos 61, sin 61) = C/2/-i(cos(0^ - 9)) = 2^cos(2j - l)(6i - cj)^). 

The rest of the proof will follow as in the case k = 21. □ 
Corollary 2.5. For k >0 and < i, j < k, 

^ 2m 

(2.11) - — — ;7fc(0i,;cos6'i,fc,sin6'i,fc)f7fc((?!)i,;cos6'j,fc,sin6'-,^fe) = {k + l)Sij. 

Z7TL ~p 1 

Proof. This follows from setting 9 = 9i^k, x = cos9j^k and y = sin0j_fc in the identity 
H2.7(l . and using the fact that 

C/fc(6'i,fc;cos6'j,fc,sin6'j,fc) = Uk{cos9i-j^k) ^ {k + l)Sij, 

where in the last step we have used the fact that C/fe(0) = (fc + 1). □ 
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Remark 2.1. The identity (|2.7|l is established for the case n = 2m. It should be 
pointed out that the result does not hold for n — 2m — 1. In fact, if n = 2m — 1, 
then following the above proof leads to an identity analogous to (|2.7|l only for 
< < m — 1 . More precisely, we have 

^ 2m- 1 

Y~ X] Uk{(l)i,; cose, sin e)Uk{(f>u; X, y) ^ Uk{9;x,y), x,y£B'^, 

v=0 



holds for < ^ < ^ 1- Indeed, in the case m = 2p it can be proved that 

2m-l 

2m 



^ 2m-l 

:;— Ura{(l)u; COS 9, sin 9)Urn{4'i^; COS (j), sin (j)) = L/„i(cos(6l-(?!))) + 2T„(cos(6'-0)), 
im ^ — ' 

where T„i{x) — cos 771-0, x = costp, is the Chebyshev polynomial of the first kind. 

There is yet another remarkable identity for Uk{- \ x, y), proved in which takes 
the summation over 9j^k instead of taking over 0^: 

Lemma 2.6. For fc > 0, 

k 

(2.12) ^Uk{9j^k;x,y)Uk{9jX,cos(j),sin(p) = {k + l)Uk{(p;x,y). 

j=o 

This identity is a consequence of the compact formula for the reproducing kernel 
of Vn{B^) in [7]- It will also play an important role in our development below. 

The two remarkable identities (|2.7I) and (|2.12|) hold the key for our new algo- 
rithms. The orthogonal polynomials in V„(-B^) are usually studied as a special case 
of the orthogonal polynomials with respect to the weight functions Wfj,{x,y) := 
(1 — x^ — y^)^~^^'^, fJ'> 0- Most of their properties are shared by orthogonal fam- 
ilies associated with for all /i > 0, and furthermore, many properties can be 
extended to higher dimensions (see |2|). The orthogonality in H2.11|l . however, is 
very special; it is not shared by any other orthogonal families associated with IV^ 
for 7^ 1/2 and there is no direct extension to higher dimensions; see, for example, 
the discussion in 



3. Fourier Orthogonal expansions 

3.1. Orthogonal expansions on the disk. The standard Hilbert space theory 
shows that any function in L^{B^) can be expanded as a Fourier orthogonal series 
in terms of V„ . More precisely, 

(3.1) L\B^) = Y,^Vk{B^): / = £proj,/, 

k=l k=l 

where proj^. / is the orthogonal projection of / onto the subspace Vk{B'^). Our 
reconstruction algorithm will be based on this Fourier orthogonal expansion. It 
is well known that proj„ / can be written as an integral operator in terms of the 
reproducing kernel of Vk{B^) in L^{B^). A compact formula of the reproducing 
kernel is given in [7]. For our purpose, we seek a formula for proj„ / that will relate 
it to the Radon transform. 
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In terms of the special orthonormal basis {Pfe : fc > 0} given in Proposition 12. 31 
we can write projj, / as 

(3.2) pioi^f^y2fj,kUk{0j,k;-), fj,k^- f{x,y)Uk{dj,k;x,y)dxdy. 

The remarkable identity (|2.7|l allows us to express the Fourier coefficients in terms 
of Radon projections: 

Proposition 3.1. Let m be a nonnegative integer. For < fc < 2m, the Fourier 
coefficient fj,k satisfies 

(3.3) f,,k = — nY, - / n^,Af;t)Uk(t)dtUk{cos{dj,k - 0.)). 

Proof. Let us denote the right hand side of H3.3() by gj^k- By the equation H2.3[l . 

1 2m 

9j,k^-^ TtXI" / f[x,y)Uk{,(j)u;x,y)dxdyUk{cos{ej^k-(l>u)) 



i/=0 



B 

2m 



/ fix^y) 7\ — j- Vc/fe(0i.;a;,y)f7fe(cos(6'j-fc - 9^^)) 

TT |^2m + l^ 



dx dy 



f{x, y)Uk{Oj,k;x, y)dx dy 

by the identity (|2.7|) . Hence, by the definition of fj^k: gj,k = fj.k- O 

A further application of the identity (|2.12l) leads us to the expression of the 
Fourier projection operator in terms of Radon projections: 

Theorem 3.2. For m > and k < 2m, the operator projj, / can be written as 
(3.4) projk f{x,y) ^ -^y]- [ TZ^^{f;t)Ukit)dt{k + l)Uk{(pu;x,y). 
Proof By (jX^ and the identity ^^J^ in Lemma [TBI 

2m 



projfc fix, y) T Y.- (•^' ^)Uk{t)dt 

fe 

X ^ C/fe(cos(6'j^fc - (j)^))Uk{OjX,x,y) 

3=0 

^ 2in ^ „1 

— V - / 7^^„(/;^)c/fe(^)d^(fc + l)^7(0,;x,y), 



2™ , _ 



which is what we need to prove. □ 
We denote the n-th partial sum of the expansion H3.1|l by Snf', that is, 



Snf{x, y) ^Y P^°jfc 



fc=0 
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The operator Sn is a projection operator from L^{B^) onto II^j. An immediate 
consequence of the Theorem 13. 21 is the foUowing result that will play an essential 
role in dering the new algorithm: 



Corollary 3.3. For m > 0, the partial sum operator S2mf can be written as 

(3.5) 52m(/; x,y)^Y.z / ^ ^-^^ (/; ^)^-(^; ^' y)di 



TT 

i/=0 



where 



^ 2m 

(3.6) ^At;x,y) ^ ^—-J2{k + l)Ukm,{(P,;x,y). 

Zm + 1 ^ — ' 

k=0 

The identity H3.5|l shows that S2mf can be expressed in terms of Radon projec- 
tions in 2m + 1 directions. This result, previously unnoticed, holds the key for our 
new algorithm given in Section 4. It should be pointed out (recall Remark l2.1|l that 
such an identity does not hold for S2m.-if- 

3.2. Summability of orthogonal expansions. Let L^'(J5^) denote the usual 
space on with norm || • ||p for 1 < p < oo and identify it with C{B^) of con- 
tinuous functions with uniform norm for p — oo. If / G U'{B'^), the error of best 
approximation by polynomials of degree at most n is defined by 

(3.7) i?„(/)p:==inf{||/-P||p:Pen2}. 

It is well-known that £'„(/)p — > for / e LP{B^) as n oo. The partial sum Snf 
of the orthogonal expansion is the best approximation to / in the norm; that is, 

\\f-S,J\\2^EM)2, fehHB^). 

However, the partial sum Snf does not converge to / pointwisely if / is merely 
continuous; see |2I and the discussion below in Section 5. 

To study the convergence of our orthogonal expansions we will introduce some 
summability methods. Such a method takes the form 

'^Oj^nSjU), a.j.n e R and ^ = 1. 

Many summability methods, for example, the Poisson means and the Cesaro means, 
have better convergence behavior (see '^). For our purpose we will use methods for 
which the operators are polynomials and we would still want to retain the property 
that polynomials up to certain degree are preserved. 

It turns out that this can be done quite easilly using a multiplier function. 

Definition 3.4. Let r be a positive integer, and let 77 e C""[0,oo). Then rj is called 
a multiplier function if 

'7(i) = 1, < i < 1, and supp?] C [0, 2]. 

Let ?7 be a multiplier function. We define an operator by 

2m / /j. \ 

S2m{f\^^v) Projfc/(a;,y). 

fc=0 
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Such an operator was used in the literature for approximation by spherical polyno- 
mials on the unit sphere, and it was been used for various other orthogonal expan- 
sions in p,U', including the expansions on the unit disk. The following theorem is 
essentially contained in 10 . 

Proposition 3.5. Let rj e C'^[0,oo) be a multiplier function. Let f G LP{B'^), 
1 < p < oo. Then 

(1) S2^f e and S^„,P ^PforPe n^; 

(2) For m > there is a constant c such that 

WSlJWp < c\\f\\p and 11/ - S'^^fWp < cE„,{f),. 

This means that the operator S'j,^ is very well behaved: it preserves polynomials 
of degree up to m and it approximates / as accurate as any polynomial of degree 
at most TO. Using Theorem 13. 21 we also have the following: 

Theorem 3.6. For to > 0, the operator can be written as 
where 

1 ^™ / fc \ 

(3.9) nit;x,y) = tE'? - (.k + l)Ukit,)Uki<P,;x,y). 

Remark 3.1. We can also use other summability methods, not necessarily prescribed 
by the multiplier function. The essence of Proposition 13.51 however, should be 
preserved. 

3.3. Fourier orthogonal expansion on a cylinder domain. Let i > and let 

Bl he the cylinder region 

Bl := ^2 X [0, L] = {{x, y, z) : {x, y) £ B\ z £ [0, L]}. 

Using the result in the Subsection 13 . II we can also get an expression for the partial 
sum operator on Bl, which will lead us to a 3D reconstruction algorithm. We 
consider orthogonal polynomials with respect to the inner product 

(3.10) {f,9)BL = -[ f{x,y,z)g{x,y,z)WL{z)dxdydz, 

where Wl is a nonnegative function defined on [0, L] with all its moments on [0, L] 
assumed finite and normalized so that J^^ WL{z)dz = 1. 

Let denote the space of polynomials of total degree at most n in three vari- 
ables. Let Vn{BL) denote the subspace of orthogonal polynomials of degree n on 
Bl with respect to the inner product H3.10() : that is, P £ VniBh) if {P,Q)bl — 
for all polynomial Q £ n^_i- 

Let pk be the orthonormal polynomials with respect to Wl on [0,L]. Let 
Uk(6j,k\x,y) be defined as in the previous subsection. 

Proposition 3.7. An orthonormal basis for Vn{BL) is given by 

P„ :== {Pn.k.j : 0< j < fc < n}, Pn,k,j{x,y, z) ^ pn-k{z)Uk{9;j.k]x,y). 
In particular, the set {P; : Q <l < is an orthonormal basis for 11^ . 



RECONSTRUCTION OF IMAGES FROM RADON PROJECTIONS 



11 



This is an easy consequence of the fact that Bl is a product of and [0, L]. For 
/ G L^{Bl), the Fourier coefficients of / with respect to the orthonormal system 
{P; : / > 0} are given by 

fi,k,j = - f{x,y,z)Pi^k,j{x,y,z)dxdyWL{z)dz, 0<j<k<l. 

^ J Br. 

Let Snf denote the Fourier partial sum operator, 

n I k 

Snf{x,y,z) =^^^fi^k,3Pi,k,j{x,y,z). 

1=0 k=0 j=0 

Just hke its counterpart in two variables, this is a projection operator and is inde- 
pendent of the particular choice of the bases of Vn(i?L). 

We retain the notation Ti^{g; t) for the Radon projection of a function g : B^ ^ 
R. For a fixed z in [0, L], we define 



(3.11) 7^0(/(•,•,z);^) := / f{x,y,z)dxdy, 

Ji{<p,t) 

which is the Radon projection of / in a disk that is perpendicular to the 2;-axis. 
The following is an analogue of Theorem l3.3l for the cylinder B^. 

Theorem 3.8. For m> 0, 

(3.12) S2mJ{x,y,z) 

= - — —T^- / n^Af{-r,w);t)<S>^{w,t;x,y,z)WL(w)dwdt 

where 



2m-k 



(3.13) <P,{w,t;x,y,z)=J2{k + l)Uk{t)Uk{^,;x,y) pi{w)pi{z). 

k=Q 1=0 

Proof. Using Proposition 13 . II and the product nature of the region, 

fi,k,j - 2 +i ^~ / T^^Mi-k^mkmtUkicosiej.k - (/-.)), 



1 --^ 



where 

fi{x,y)^-l f{x,y,w)pi{w)WLiw)d'w, l>0. 
Jo 

Substituting this expression of fi^kj into the formula of S2mf and using the identity 
in Lemma EEl we then obtain 

S2mf{x,y,z) = —-—Y,- / 

v=0 -^^^ 1=0 k=0 

X {k + l)Uk{(f>u; x,y)pi_k{z) 

2m ^ ^1 2m 2m— k 

= r^Y.- T^<t>Afi:t)pi{z){k + l)Uk{c^,-x,y)Uk{t)dt, 

v=o ■^-'-k=o 1=0 
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where in the second equahty we have exchanged the two inner summations. From 
and the definition of /; it is easy to see that 

= - / n^,Af{-r,w)-t)pi{w)WL{w)dw, 

Jo 

from which the proof follows upon rearranging terms in the summation. □ 

Remark 3.2. Clearly one can also consider summability of orthogonal expansions 
on Bl. For example, one can define the operator with multiplier factors just as in 
the case of orthogonal expansion on B^. Wc shall not elaborate. 

4. New Reconstruction algorithms 

4.1. Reconstruction Algorithm for 2D images. The identity (|3.5|l expresses 
S2mf in terms of the Radon projections TZr/,^ (/; t) of 2m + 1 equally spaced angles 
(pi,, < 1^ < 2m, along the circumference of the disk. These Radon projections 
are defined for all parameters t. In order to make use of the Radon data from the 
parallel geometry, we will use a quadrature rule to get a discrete approximation of 
the integrals 

n^M\t)Uk{t)dt = j'^ '^d^Uk{t)VT^dt 

in (|3.5|l . The result will be our algorithm. 

If / is a polynomial then the equation (|2.5(l shows that Tl^{f\ t) /Vl — t"^ is also 
a polynomial. Hence, we choose a quadrature rule for the integral with respect to 
\/l — on [—1, 1]. Let us denote such a quadrature rule by Ing'-, then 

(4.1) - / g{tyi-t^dt « ^^9{tj) In(ff), 

where ti, . . . ,t„ are distinct points in (—1, 1) and Xj are real numbers such that 
Sj=i — 1- If equality holds in Ij4.1(l whenever g is a polynomial of degree at 
most p, then the quadrature rule is said to have the degree of exactness p. 

Among all quadrature rules with a fixed number of nodes, the Gaussian quadra- 
ture has the highest degree of exactness. It is given by 

(4.2) i git)VT^'dt = ^ i: sin^ -^.g f cos := I^ig) 

for all polynomials g of degree at most 2n — 1; that is, its degree of exactness is 
2n — 1. Note that j7r/(n + 1) are zeros of the Chebyshev polynomial C/„. 

Using quadrature formula in (|3.5|) gives our reconstruction algorithm, which 
produces a polynomial A2mf defined below. 

Algorithm 4.1. Let the quadrature rule be given by H4.1|l . Form > and {x,y) G 
B^ 

2rn n 

(4-3) A2m.{f;x,y) = ^^7^0„(/;^J)^J,^(a;,y), 

where ^ 

Tj,v{x,y) = ^ <^„{tj;x,y). 

2(2m+ 1)^1 -t2 
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For a given /, the approximation process A2mf uses the Radon data 

{n^Af;tj):0<iy<2m, 1 < j < n} 

of /. The data consists of Radon projections on 2m + 1 equally spaced directions 
along the circumference of the disk (specified by (f>^) and there are n parallel lines 
(specified by tj) in each direction. If these parallel Radon projections are taken 
from an image /, then the algorithm produces a polynomial A2mf which gives an 
approximation to the original image. 

The polynomial A2m is particularly handy for numerical implementation, since 
one could save T^^^ in a table beforehand. This provides a very simple algorithm: 
given the Radon data, one only has to perform addition and multiplication to 
evaluate A2m in (|4.3(l to get a reconstruction of image. 

A good choice of the quadrature rule is Gaussian quadrature. We choose in 
particular n = 2m so that the nodes of the quadrature rule H4.2|l becomes tj = 
cos9j^2m = cosj7r/(2m + 1). In this case, our algorithm takes a particular simple 
form. 

Algorithm 4.2. For to > 0, {x,y) e B^, 

2m 2m 

(4.4) A2m{f;x,y) = ^^7^0J/;cos6'J,2m)^J>(x,?/), 

i/=Oi=l 

where 

^ 2m 

(4.5) Tj^,{x,y) = — — ^2 ^(fc + l)sin((fc + l)0,,2„OC/fc(0.;x,?/). 

' k=0 

Remark 4.1. It is tempting to choose n — 2m + 1 in the Gaussian quadrature. 
In fact, such a choice has one advantage: the operator A2m will be a projection 
operator onto IIjj^. We choose n = 2m so that (pi, — 2m+i ^-"^"^ ^i,2m — 2m+i have 
common denominator. It also turns out that this choice works perfectly with the 
fan beam geometry of the projections, which will be reported elsewhere. 

The convergence of the algorithm will be discussed in the following section. Here 
we state one result that is a simple consequence of the definition. 

Theorem 4.3. The operator A2m in Algorithm \4. l\ preserves polynomials of degree 
a. More precisely, A2mif) = f whenever f is a polynomial of degree at most a. 
In particular, the operator A2m in Algorithm \4. '(^ preserves polynomials of degree at 
most 2m — 1. 

Proof. The polynomial A2mf is obtained by applying quadrature rule H4.1|l to Ij3.5(l . 
If / is a polynomial of degree at most a, then so is TZ^^^ (/; t) / \/l ~ t'^ for every v. 
Furthermore, the polynomial ^,y(t;x,y) is a polynomial of degree 2m in t. Hence, 
when we apply the quadrature rule H4.1|l to the integral l|3.5|l . the result is exact. 
Therefore, A2mf = S2mf - /• □ 

If the quadrature rule is the Gaussian quadrature in H4.2() . then the polynomials 
preserved by the operator A2m have the highest degrees among all quadrature rules 
that use the same number of nodes. Such a choice will ensure better approximation 
behavior of A2m- 
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Remark 4.2. Using the angles (jfiy, another projection operator, call it 1/2™/, has 
been constructed in ^ based on the parallel geometry. For almost all choices of 
{tj S (—1, 1) : < J < m}, the operator is the unique polynomial of degree 
2m determined by the conditions 

7^0^ (J2m.f ; U) = 7^^^ (/; U), 0<iy< 2m, < j < m. 

However, the construction of requires solving a family of linear system of 
equations whose coefficient matrices, depending on the choice of tj, appear to be 
badly ill-conditioned. 

4.2. Reconstruction algorithm for 2D images with a multiplier function. 

Instead of using (|3.5(l . we can also start from the summabilty with the multiplier 
factor H^18|l . The result is another reconstruction algorithm. Here we state the 
resulting algorithm only for the Gaussian quadrature H4.2() with n = 2m. 

Algorithm 4.4. For to > 0, {x,y) G B^, 

2m 2m 

(4.6) A^2mif;^,y) = EE^*^(/;co«^j-.2m)T;';.(x,2/), 

where 

1 / fc \ 

" (2^+1)2 E ^ [-) + 1) ^'^((^ + l)^i.2„OC//c(0.; X, y). 

For a given /, the approximation process A'lmf uses the same Radon data of / 
as A2mf ■ It also has the same simple structure for numerical implementation. Its 
approximation behavior appears to be better than that of A2mf- We conclude this 
subsection with the following analogous of Theorem l4.3l 

Theorem 4.5. The operator preserves polynomials 0/ degree to. More pre- 
cisely, A'lmif) — f whenever f is a polynomial of degree at most m. 

We can obtain other reconstruction algorithms using different summability meth- 
ods; see, however. Remark l3. II 

4.3. Reconstruction algorithm for 3D images. In order to get a reconstruction 
algorithm for 3D images on the cylinder region Bl, wc choose the weight function 
to be 

Wl{z)^-^=^=, ze[0,L], 

which is the Chebyshev weight function on the interval [0, i], normalized so that 
its integral over [0, L] is 1. Let be the Chebyshev polynomial of the first kind. 
Define fu by 

fo(z) = 1, ffc(z) = V2Tk{2z/L - 1), k> 1. 

The polynomials T^, are orthonormal polynomials with respect to on [0,L]. 

To obtain an algorithm using parallel Radon projections, we start from (|3.12|) 
and apply quadrature rules on its integrals. For the integral in z, we use the 
Gaussian quadrature for Wl- Set 

{2i ~\- 1)tt l + cosfi„ 
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where Zi are zeros of Tn{z); then the Gaussian quadrature on [0, L] takes the form, 



(4.7) - f\{z) . , = -Eg(^0 



Jo y/z{L - Z) 



1=0 



which holds whenever g is a polynomial of degree at most 2ri — 1. For the integral 
in t, we could use the quadrature rule H4.1|l . For simplicity, however, we will only 
use the Gaussian quadrture (|4.1|l . 

This way we get an algorithm for reconstruction of images on ■ The algorithm 
produces a polynomial B2m of three variables as follows: 

Algorithm 4.6. Let — TZ^^{f{- ,-, Zi); cos 9 j_2m) ■ Form>0, 

2m 2m n— 1 

B2nif{x,y,z) := ^^^j^j^,T^j^,{x,y,z), 



wh- 



ere 



v=0 j=l 1=0 

1 



Tv,j,i{x,y,z) = — — -—<^^{z,,cos6j^2m]x,y,z). 

n[2m +1) 

For a given function /, the approximation process B2m uses the Radon data 

{7^0„(/(•,•,z,;cos6lJ■2„O : < 1/ < 2to, 1 < j < 2to, < ^ < n - 1} 

of /. The data consists of Radon projections on n disks that are perpendicular 
to the z-axis (specified by Zi); on each disk the Radon projections are taken in 
2m + 1 equally spaced directions along the circumference of the disk (specified by 
(/)y) and 2m parallel lines (specified by cosipj) in each direction. We can use this 
approximation for the reconstruction of the 3D images from parallel X-ray data. In 
practice, the integer n of z direction should be chosen so that the resolution in the 
z direction is comparable to the resolution on each disk to achieve isotropic result. 
The following theorem is an analogous of Theorem 14. HI for Bl- 

Theorem 4.7. Ifn> 2m, then the operator B2m in Algorithm \4-(A preserves poly- 
nomials of degree 2m — 1. More precisely, B2mif) = / whenever f is a polynomial 
of degree at most 2m — 1 . 



Proof. If / is a polynomial of degree at most 2m — 1, then 7?.0(/(-, •, w); t)/\/l — t"^ 
is a polynomial of degree at most 2m — 1 both in the t variable and in the w 
variable. The function $y (w, t; x, y, z) is of degree 2m in both the t variable and 
the w variable. Hence, when we use the quadrature rules for S2mf in 13.12|l . the 
result is exact if the quadrature rules are exact for polynomials of degree 4m — 1. 
For the quadrature rule 1)4. 7|l this holds if n > 2m. □ 

Remark 4.3. In the z direction, we choose the weight function (z(l — z))~^/'^ in- 
stead of constant weight function. The reason lies in the fact that the Chebyshev 
polynomials of the first kind are simple to work with and the corresponding Gauss- 
ian quadrature 1)4. 7|l is explicit. If we were to use constant weight functions, we 
would have to work with Legendral polynomials, whose zeros (the nodes of Gaussian 
quadrature) can be given only numerically. 
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5. Convergence of the reconstruction algorithm 

In this section we study the convergence behavior of the reconstruction algorithm 
for 2D images and we work with A2m in H4.4|l . for which the quadrature rule is 
chosen to be Gaussian quadrature. 

5.1. Convergence of 2D reconstruction algorithm. Let us consider the uni- 
form norm || • ||oo of continuous functions on B^. Convergence in the uniform norm 
guarantees the pointwise convergence of the reconstruction. First we give a formula 
for the operator norm. Let us denote by 

i^j Oj^2nr = j7r/(2m + 1), < j < 2m. 

Proposition 5.1. Let ||.42m||oo denote the operator norm of A2m o-s an operator 
from C{B^) into C{B^). Then 

2m 2m 

(5.1) ||^2m|loo = max Ajn{x,y), A™(x, y) := sin?/>j|Tj.,,(x, ?/)|. 

("^^^^^ .=o,=i 

Proof. By the definion of TZ^{f;t) we evidently have 

< Vl-t^f\\oo, O<0<27r, -l<t<l, 
as seen from the second equality of (|2.2|) . It follows immediately that 

(5.2) |^2m(/;a;,y)| < ||/||ooA™(x,?;), where {x,y)eB^. 

Taking maximum in both side proves that ||yl2m||oo < ''^^^(x,y)eB^ \-^m{x,y)\. To 
show that the equality holds, let (xq, j/o) be a point in at which A^ix, y) attains 
its maximum over B^ . Recall that l{0,t) denote a line segment (|2.1|l inside B^. 
Let E denote the set of intersection points of any two line segments I{4>i/ , cos ) 
and I{4>fi, cosipi), 

S := {{x,y) : /(0^, cos 7/)^) n /((/)^, cos Vi), 7^ 

The set contains only finitely many points. Let e > be small enough so that a disk 
centered at a point in E of radius e contains no other points in E. Let denote 
the union of all such e disks. We construct a function e C{B^) as follows: 

f{x,y) = sintljj signTj,i,{xo,yo), {x,y) £ I {(j>^, cos (t)j) \ (/((/)^, cos 0^) n E) 

foraU j,j/and H/eHoo = 1- Then 7?.^„(/£, cos -0^) = sign Tj>(xo, yo) + Cj>(r for some 
costant Cj.i/. Since there are only finitely many points in E, this shows that 

||-42m||oo > |-A2m/e(2;o,yo)l = Am(a;o,yo) - ce = max A,n{x,y)-ce, 

where c is a constant that depends on m. Taking e ^ completes the proof. □ 

In the following we shall use the notation A B to mean that there exist two 
positive constant ci and C2 such that ciA < B < C2A. 

Recall that A2mf is obtained from the partial sum S2mf of the Fourier orthogo- 
nal expansion, upon using Gaussian quadrature. It is proved in j^j that the operator 
norm for the partial sum operator S2m is 

(5.3) ||S'2„i||oo ~ "I- 

One would expect that the operator norm of A2m is worse than 0{m) due to the 
additional step of using Gaussian quadrature. Our main result in this section shows 
that the norm of A2m does not grow much worse. 
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Theorem 5.2. For A2m defined in (|4.4|l . 

||-42,„||oo ~ mlog(m + 1). 

This shows that the price we pay for using Gaussian quadrature to get to .42m 
is just a log(TO+ 1) factor. This theorem wiU be proved in the foUowing subsection. 

For / e C{B^), the quantity £'ri(/)oo defined in (|3.7|l denotes the error of best 
approximation of / by polynomials of degree at most n. It is proved in |10| that if 

/ e C^'-iB^), r e N, then 

(5.4) S„(/) < cn-2'-||2?''/l|oo, n > 0, 

where 2? is a second order partial differential operator defined in (|2.4|l . 
As a consequence of Theorem 15. 21 we have the following result: 

Theorem 5.3. If f & G^{B^), then A2mf converges to f uniformly. In fact, let r 
be a positive integer: then for f G C^^{B^), 

l|A„/-/||oo<ci^^5%±ll||PVI|oo. 

Proof. Let p be the best approximation for / from n^^. By the definition of the 
operator norm, the fact that .42m preserves polynomials of degree up to 2m — 1, 
and the triangle inequality we see that 

M2m/-/||oo < \\CA2,n{f-p)\\oo + \\f-p\\oo 

2m||oo)-£'2m— l(/)oo !i Cm log(m + l)-£'2m— i (/)oo 

from which the stated inequality follows from (|5.4I) . □ 

This theorem shows that the Algorithm l4.2l does converge whenever / is a C^{B^) 
function. In other words, if the original image is smooth then the reconstruction 
algorithm 14. 21 converges to the image pointwisely and uniformly. The speed of the 
convergence depends on the smoothness of the function. 

The algorithm with multiplier function will likely have better convergence be- 
havior (recall Proposition 13. 5|l . but the estimate is more difficult to establish. We 
will report results along this line in future communications. 

5.2. Proof of Theorem 15.21 lower bound. First we need a compact formula for 
the functions Tj in (|4.5|) . The notation 

cos Ov{x,y) = X cos <p^ + y sin (fii,, (p^, = ^nvj (2m + 1). 

will be used throughout the rest of this paper. 

Proposition 5.4. For Q <v < 2m and I < j < 2m, 

(5.5) (2m + irra-, y) = ' ""^i^^ " ^-/^^"^^"^^('^"^y^' ^^^^ 

2{cos i,{x, y) — COS ipjY 

- ( lV::mi'j (^"^ + ^)^2m (cos 6l^(a;, y)) 
[ ) sm^fjj 2(cos6'^(x,y) - cosV-j) 

Proof. To derive the formula, we start with the elementary trigonometric identity: 

V^ri.^i^...^^i.^i^m - -1 + (2m + 2) cos((2m + 1)9) - (2m + 1) cos((2m + 2)9) 
2J^k+l)cos((k+l)&) - 4sin2(0/2) ■ 
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Let 6 = 9u{x,y) in this proof. We apply the above identity to 

2m 

sin6l(2m + lfTj^^{x, y) = ^(fc + 1) sin((fc + l)V'j) sin((/c + 1)61) 

fe=0 

= 2 + 1) + - ^J-)) - cos((fc + 1)(^ + V,))] 

fe=0 

and combine the two terms together as one fraction. The denominator of the 
fraction is 

2 . 4sin2 sin^ = 2(cose - cos^,-)' 

and the numerator of the fraction is 



N{e) := sin^ t^hj,m{0 - i^j) - sin^ t^hj,m{0 + V',) 

where 

hj,m{d) = -1 + (2m + 2) cos((2m + 1)9) - (2m + 1) cos((2m + 2)6*). 
We write N{9) as a sum of two terms, 

Ni9) = Ni{9) - (2m + 1)N2{9). 
For the first term wc use the identities 

cos((2m + 1)(6' ± V'j)) = (-1)^' cos((2m + 1)6'), 
and 2sin^(6'/2) = 1 - cos 6* to get 

Ni{9) :=sin2 ^-^(-1 + (2m + 2) cos[(2m + 1)(^ - ipj)]) 

- sin^ ^~^^ (-l + (2m + 2) cos[(2m + l)(6l + ipj)]) 

=(-1 + (2m + 2)(-l)J' cos[(2m + l)6'])(cos(6' - tpj) - cos(6' + ^j))/2 
={-l + {2m + 2){-l)^ cos[{2m+l)9])sm9simjjj). 

For the second term, we use 

cos((2m + 2)(6' ± ^Jj)) = (-1)^' cos((2m + 2)6* ± V-j), 

2sin^(^?/2) = 1 — cos^ and the addition formula of the cosine function to obtain 

N2i9) := sin^ cos((2m + 2)(6l - Vj)) - sin^ ^ ~ cos((2m + 2)(6l + Vj)) 

- ^ ' [cos((2m + 2)61 - Vj) - cos((2m + 2)6* + tpj)] 



2 
2 



+ ^7^[cos(6' - Ipj) cos((2m + 2)9 - -0^) - cos(6' + V'j) cos((2m + 2)6* + ipj) 



= (-1)^ sin V'j sin((2m + 1)6*) + ^^^^h cos(6' - ipj) sin(6» + ipi) sin((2m + l)6i) 

+ cos(6i + Ipj) sm{e - ipi) sin((2m + 1)^)] 
= (-1)-'' sin'^j sin((2m + 2)61) + (-1)-''+^ sin((2m + 1)6*) sin V'j cos V'j, 
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where we have used the double angle formula for sine in the last step. Using the 
addition formula sin((2rn + 2)6) = sin((2TO + 1)6) cos6 + cos((2to + l)^) sin0, we 
obtain 

N2{9) ^ (-1)^ smil;j[sm{{2m + l)6){cos6 ~ costpj) + cos{{2m + 1)6) sin 9]. 
Putting the two terms together we obtain 

N{6) = - sin6'sin'(/'j + (-!)■'' sin V'j [(2m + 2) sin6'cos((2m + 1){6)) 

- (2m + l)(sin((2m + l)6')(cos6' - cosipj) + cos((2to + 1)6') sin 6*)] 
= - sin6lsinV'j[l - (-1)^' cos((2m + 1)6*)] 

— (2m + 1)(— 1)^ sin^j sin((2m + l)^)(cos6' — cos^'j)- 
Consequently, we have proved that 

-sine'sinV'j[l - (-1)-'' cos((2to + 1)6*)] 



sin6'(2m + l)^rj>(a;,y) 



2(cos6' — cos ■0^)2 
(2m + 1)(-1)^' sinipj sin((2m + 1)6*) 



2 (cos 6* — cos-0j) 

Hence, using the fact that T„(cos0) = cos n6 and Un{cos6) — sin(n+ 1)6/ sm6, we 
conclude that 

(2m+l)2T,-,(a;,y) 

_ -smiPjjl - (-l)Jr2^+i(cos6>)] _ (2m + 1)(-1)^' sinV>j;72^(cos6>) 
2(cos6' — cosi/jj)^ 2(cos6' — cos-i/ij) ' 

which completes the proof. □ 

Throughout the rest of this section and in the following section, we will use the 
convention that c denotes a generic constant, independent of / and m, its value 
may change from line to line. The elementary facts 

-t<smt<t, forO<t<f, and cos a - cos/3 = 2 sin j sin j 

will be used repeatedly without further mention. 

We now use the expression (I5.5|l to derive the lower bound for the estimate. 

Proposition 5.5. 

/ TT TT \ 

M2m||oo > Am I COS ^ ^ , sin 2 ) - cmlog(m + 1). 

Proof. Let x = cos and y = sin^;;^. Then cos6^{x,y) = cos 2m+f ' 

that sin(2m + l)6i,{x,y) — 1 and cos(2m + l)6i,{x,y) = 0. Let 0^ ~ ''^2m+i^^ - 
Then by JCTl . 



(2m+l)2r,,,(x,y) = -^ 



2m + 1 



{cos 6^ — cos -0^)2 sin (cos — cos^pj) 

We will use the fact that < 6,^ < tt and < tpj < 7r/2 for Q < v,j < m. 
Furthermore, for < ^i, < 7r/2 and < ■0 < ''i'/2 we have 

sin ibi 2 sin cos % ibi 

rx5- = 1 < 2 cos < 2, 

sm — sm ■'-^ — ^ 
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which also holds for -n 12 <Qy < tt and < ip < 7t/2, since then7r/4 < ^^^^ < 37r/4 

2 



and sin ^^'^^•^ > ^/2/2. Hence, we have 



1 x ^ x ^ Sm %f}j 1 X - X ^ 1 

(2m + 1)2 2^2^ (cos 6*^ - cosV-iO^ ~ (2m + 1)2 ^ ^ cin2 



< 



EE (2. 1/2)2 ^^"^- 



Consequently, it follows that 

^ mm 
Am{x,y) > 7^——T^Y.J2''''^^l\^oAX:y)\ 

, mm • 2 / 

^ 1 ■^-^ ■^-^ sm t/Jj 



tEE' 



2m + 1 ^ ^ I sin (cos 6*1/ — cos'!/',)| 
The sum in the last expression is bounded below by 



8 1 V^l 

^2 2^+1^^^0,10,-^,1(0, 



EE 



2 



^^1 - 1/2) |2j^ - J - 1/2| {21. + J - 1/2) 

^i:^^^^^ 1 _ 1 - 1_ 1 :^^m/2-, + l 



7]-3 2i^ — 7 TT^ -^-^ -^-^ 7 TT^ -^-^ 7 

j=v •' v=l j=l ■' j = l ■' 

m/2 

1 /m \ X ^ 1 m , , , , 

^^(y + i)E--^^^-i°g(- + i)- 

This completes the proof. □ 



5.3. Proof of Theorem l5.2l upper bound. We will also use the expression (|5.5f) 
to estimate Am{x,y) in (|5.1|) from above for all {x,y) e B^. In the following we 
write A{x,y) = A™(x,?/). 

It is easy to see that the function A(x, y) is invariant under the dihedral group 
^2m+i; that is, it is invariant under the rotation of a angle for v = 0,1, . . . , 2m. 
Hence, it suffices if we establish the estimate assuming (x, y) is in the wedge 

7r/2 

:= {ix,y) : X = r cos = r sin 0, r > 0, |0| < £m = 7; — • 

Zm + 1 

Note that the set Tm is symmetric with respect to the y axis. 

We start with a number of reductions. The fact that (j)2m+i-i/ — 2Tr — c/)^ shows 
cos02m+i-i/ = COS0, and sin (/)2m+i-iy = — sin^,, which implies that 02m+i-i/(a;, y) = 
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9v{x,—y). Hence, using that Oo{x,y) — x, we can write 



2m 



|ro,.(2:,y)| +^ (|r,-.(a:,y)| + -y)|) 



!/=l 



u=l 

Since €; Fm, we only need to consider the sum over H^{x,y). 

The equation H5.5|l shows that Tj^^, is naturally split as a sum of two functions, 



(5.6) 
(5.7) 



(1) -sin^j[l - {-iyT2m+i{cos9^{x,y))] 



T. 



(2) 



2{cos 6i,{x,y) — cosipj)^ 



J. I' 



-(— l)-* sin 



{2m + l)U2rn{cose^{x,y)) 
2{cos9y{x,y) — coa-ipj) 



We shall denote the corresponding splitting of H^{x, y) by HI, {x, y) and the split- 
ting of A(a;, y) (only the sum over H^{x, y) by A''*-'(a;, y). Thus, 

m 2m 

AW(x,y) :=^i7«(x,y) and i7«(x, y) ^ sin^,|i;(;3(x, y)|, z = l,2. 

Next we split H^,^^ as two sums; the first one is over j = 1, 2, . . . , m and the 
second one is over j = m + 1, m + 2, . . . , 2m. Since ip2m+i-j = tt — tpj , we have 
cos V'2m+i-j ~ — cos ijjj so that we can write 

(5.8) Hi^H^,y) - Hl%,y) + H^},{x,y), 

where, using that 1 — (—1)-' cos(2m + 1)9 = 1 — cos(2to + 1){9 — tpj), 



(5.9) 



"■^^ (2m +1)^^ 

j=i 

rn 

"■'^ (2m +1)^^ 

1=1 



1 ~ cos{2m + 1){9 - -ipj) 



2{cos9,^{x,y) — cos -0^)2 
1 ~ cos{2m + 1){9 - ipj) 



2 (cos 0^ (a;, y) + cosipj 



and, using the fact that (—1)^ sin(2r7i + 1)6' = sin(2m + 1){9 ~ tpj)! 



(5.10) 



4a (^,y) 



Let us denote the corresponding split of A'^*^ by A^*-*; that is, 

m 

A«:=Ai')+A«, where Af z,j = l,2. 

1^=1 

We only need to estimate one of the two terms. To see this, let us define 
0^ (2i/- l)7r/(2m + 1), 1 < < m. 



2m 



- sin^ tPj 



sin(2m-|- l)(^9„(^x,y) - tpj) 



2sin6',y(a;,y)(cos6'^(a;,y) — cos-^j) 
sin(2m + l){9t,{x, y) - tpj) 



2 sin 6*1. (a;, y){cos9i,{x, y) + cosi/'j) 
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Then COS 1-1/ = cos(7r— 0^) = — cos0^ and sin0m_,y+i — sin^jy. Consequently, 

cos9rn~i^+iix, y) — — X COS + ys\i\(f)i, — — cos 6*1/(2:, — y), 

where 6^{x,y) — xcos.4>^ + ysm^i,. We will also use the notation and A^'^ 
when 9^{x, y) is replaced by 9v{x, y) in hIj\. It then follows from (|5.1U|) that 

H^m-v+iA^iy) = Hl'l{x,~y), l<v <m, 



and, consequently, K2\x,y) = K^f' {x,—y). Hence, the estimate for Aj''' will be 

7r/2 



(2) 

similar to the estimate for h\ . In fact, set 



Tm := {{x,y) : X ^ rcos(j>,y ^ rsin^, r > 0, \tt ~ (f)\ < e„}, = 



2to + r 



then the estimate of A2 (x^y) over Tm will be exactly the same as the estimate of 
A^*''(a;,y) over r,„. Thus, we only need to estimate one sum, which we choose to 
be A«. 

We use the equations H5.9|l and H5.1()(l to carry out the estimate. The inequality 

sin(n + l)t 



(5.11) |f/„(cost)| 



<n+l, < t < TT, 



sin^ 

will be used several times. The estimate is divided into several cases. 
Lemma 5.6. There exist constants ci and C2 such that 

Hq^\x,v) <ci and H^^^ (x, y) < C2 + log(m + 1) 

for {x,y) e r,n, where ci < n^{2 + 77^/12) and C2 < + 1. 

Proof. Since cos9o{x,y) = x and {x,y) e r,„, we can write 9 = 9^{x,y) with 
0<9< 7r/2. 

(2) 

We estimate the term Hi, ' first. Since cosipj > for 1 < j < to and cos 6* > 0, 
it follows from ifOTjl that 

"•2 y> -2^ cos iP, - ^ COS V,n-, ^ sin (■''+^/^)" 

TT^ 2TO+1 / 1\, , 

%E(-TT72KH"^+2ji°s(^" + ^)- 

f 2) 

For the term ( , we need to consider the position of 9 relative to that of ipj , 
which is divided into several further cases. 

(A) If < 6* < e„i ^ 7r/(4TO + 2), then by (|5lT|) 



• 2 / m 1 2 m 1 2 

-"0,ll2;,2/j . ^ . ^ +e S TT 2^ ,2 fl2 - ,/,2 =-2 

4 sm -^^^ sm -^^^ - , V, — t' - , V, — 



™ -2 / 1 ™ 1 \ 

- ^tt - 1/4) - 4 {p - 1/4) J - 
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(B) li ^pi — Em < < ipi + Em, whcre 1 < ^ < to, then using the fact that 



sm(2m+ ■(/') 



sm 



< 2 



sin(2m + l)^ 



< 2m+ 1 



we obtain that 

r(2) 



sin^ V/ 



i-i 



sin^ "0) 



<- 



8 9{^i + 9) 8(2m + 1) 



^ — 1 m 

E+E 



±2+1 
2 2 



E 



E - 



.1-1/2 8(2m + 1)9 yf^^ I- J- 1/2 ^.^^^^ 
The first sum in the last expression is bounded by 



J -I- 1/2 



, .^^ /-j-1/2 - 8{rTj2) h ^-.^--1/2 - T^^^S"^ + 



and the second sum is bounded by 

•2 m 



< 



2 m — l . , 

^ J +1 



8(2m + 1)9 ^.^^ J - ; - 1/2 - 8(/ - 1/2) f^^ j - 1/2 



E 



< 



9 — / 

TT m 



(logm + 2). 



Putting these estimates together and use (|5.8|) , we complete the proof for Hq^'^ . 

The proof for Hq^^ does not need require the splitting argument. By definition 
and H5.6|l . 

2rn 



j=i 

^ 2m 

j=i 



l~cos{{2m+l){9~iPj)) 
Ssin^^sin^^ 

sin^ ((2m + 1)^) 

^ ^sm ^ 



Let 9 be fixed and |0 — "0^1 < £m for some fc. Then by (|5.11|l . 



+ 



SVC? ipj 



sin tpk 

Asm' 2^ ' (2m+l)2^^4sii,2^gi^2fl^ 



fe-i 



< 



xh-E 



E 



2 ' ^ (j-fc- 1/2)2 



< 



This completes the proof. 



□ 
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Putting all pieces together, it is readily seen that the proof of the Theorem 15.21 
follows from the conclusion of the following lemma. 

Lemma 5.7. There exist constants ci and c2 such that 

A\^\x,y) < ci (m + 1) and A^^''(a;, y) < C2 (m + 1) log(TO + 1) 

for {x,y) e r„, where ci < Tr^/gG + 77rV24 + 1/2 and C2 < (3/2)7r2 + tt + 1. 

Proof. Again we consider the estimate for the sum h.\ first, which is more difficult 
than the estimate for A^^-*. Using the polar coordinates x = rcos(f> and y = rsincj), 
we can write cos9^{x,y) — rcos{(f>~ (j)^). For {x,y) G F^, \4>\ ^ Sm, so that 

(2:^-1/2)^ ^ (2z.+ l/2)^ 

2m + 1 - 2to + 1 ' 

from which we conclude that cos 6*1^(2;, y) > for 1 < 1/ < m/2 and cos6',y(a;, y) < 
for m/2 + I < < m. 

We shall break the sum in A^*'' into two parts, 

A«(x,y)=L«(a;,y) + L«(x,y), 

where the first one has the sum taken over 1 < v < m/2 and the second one has 
the sum taken over r7i/2 + \ < v < m. 

Case 1. We consider the estimate of L2\x,y). Note that cos9u{x,y) < and 
cos^pj > for the terms in L2\x,y). If sm6j^{x,y) > \/2/2, then 



1 


2m 4 


1 


1 




2m 4 


1 


2m 4 


1 



i/=m/2+l i=l 



sin((2m + l)6'j.(x,y)) 



2 sin 01, {x, y) {cos 6iy{x, y) — cos ipj) 



^ ^ ^ l_ sin ^ 1 ^ 1 



^ J2 cos ^Jj ~ 4V2 ^ sin ^^+^/'^'>'' 

< r=- > < — ^ m + - logfm + 1). 

- iV2 f^^j + 1/2- 4V2\ 2 J ' 

If sin 6*^(0;, y) < \/2/2, then — cos0y(a;, y) — cos(7r— 0^(a;, y)) > a/2/2. Furthermore, 
since — cos9v{x,y) — —rcos{(f)y — (j)) < — cos(0j/ — (j)) for m/2 + 1 < < m, we 
have 0„{x,y) < cj),j - (j) and | sin6'i.(x, y)| = sin(7r - 9„{x,y)) > {2/tt){tt - + <!>)■ 
Hence, for sin6'i,(x,y) < \/2/2, 

-{2)f . ^ 1 1 2sm'' Vi 



^2'^^^'^^- 2m + l ^ E^lsin^.C 

V =—y 



TT 1 1 ^ 2m + 1 

i^=m/2+l " ^ " ' ^ - V - [/=! 



2V2 , ^ ^ - 0. + 2V2 + 1 
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The case {^^ v) is again easier. We have 



^ mm 



(2m + 1) 



i/=m/2+l j=l 



sin2((2m + l)M£|)^) 



2(cos6'i,(x,y) — cosipjY" 



< 



■E 



sin ipj m 



< 



4(2m + 1)2 cos2 V'j ~ 16 ^ {j + 1/2) 



TO X ^ 1 

1 R (n J- ^ 



Case 2. The estimate of L^i\x,y) under the assumption that 
< r < cos{tpm — £to) = sin7r/(2m+ 1), 

where cos9v{x, y) = rcos{(j) — (j^v) in the polar coordinates for (x, y). The fact that 

cos6i,{x,y) < r impHes that 9^{x,y) > ipm — £m = 7r/2 — 2£m, so that sin 6*^(2;,?;) > 



1/2. Furthermore, 7r/4 < ^ 

\/2/2 > 1/2. Consequently, 

m/2 



< 37r/4, which implies that sin 



> 



lit 



u=i j=i 



sin((2TO+l)(6l^(x,y)-Vj)) 



2sin '-<"f-'^- sin '-("f+'^' 



m/2jri— 1 ■ 2 / 

1 V V^j I V 

1 ^ |g4^,y)-'A.I ^ ^ 

iy=l j=l 2 



sin^ ^pr, 



2m- 



\eA^,y)+i'm\ 
2 



< 



m — 1 



1 



+ 2m 



2 ^ m - 7 - 1/2 
Similarly, we have the estimate 

^ m/2 m 

^ ' iy=l 7 = 1 



2(2to + 1) 7I'/2 - Vi - 2£„ 

m— 1 ^ 

+ 2m < m(logm + 2). 



^tE 



sin^ ((2m + l) ^-^"f-^^ ) 



4 sm ^ 2 sm ^ '"^T^yj 



/2m-l 



< 



^/2 



_ sin^ 
4(2m + 1)2 (^.(a;, y) - V,)^ ^ ;^ 4sin2 



< 



TO 



E 



TO TO / TT 
+ ^ < 



+ 5 



2 ^ (m-i - 1/2)2 2 - 2 V 6 
Case 3. The estimate of L'f^ (x, y) under the assumption that 

r > cos{^pm — £m) = sinTr/ (2to + 1), 

where cos 6^{x,y) = rcos((/) — (pi,) in the polar coordinates for {x,y). Note that 
cos 9,^{x,y) > for the terms in L''^\x,y). In this case, cos 9,y{x,y) — cos^j can be 
zero. For a fixed r let I be the index such that 

r = cos^, \9 — tpi\ < £m, 1 < I < m. 



26 



YUAN XU 



We split the summation over j as two sums, one over 1 < j < I and the other over 
I + 1 < j < m. This leads to a split of the sums in L^\x, y), 

m/2 I m/2 rn 
v=l j=l i/=lj=/+l 

We consider the two cases separately. 

(A) The estimate for L^\. For j < I, r = cos6 < cos{ipi — e) < cosipj, so that 
cos 9v{x,y) <r< cos ipj, and consequently, 

I cos9u{x,y) — cosV'jl > cosipj — cos6v{x,y) > cosipj — cos{ipi — Sm) > 0. 

Furthermore, the fact that cos9v{x,y) = rcos(0— (pv) < cos('^; — Sm) implies 1 
9v{x,y) >ipi- Em- Therefore 



Furthermore, the fact that cos9v{x,y) = rcos(0— (pv) < cos('^; — Sm) implies that 
9v{x,y) >ipi- Em- Therefore 

7n/2 I . 2 , 

1-1^ ' - 2m + 1 4| sin^,(x, y) sin M-:^ sin 

~ 8(2to + 1) V'; - £m ^ (V'i + V'j - em)(V'z - V'j - £m) 

The estimate for is similar, 

^(1). 1 \^ sm Ipj 

- (2m + 1)2 A. 8| ^j^2 M£|)^ sin^ 

TT^m 1 TT^m /tt^ 



(B) The estimate for i^*2- I < j < m, r — cos 9 > cos(tpi + e) > cos^'j-; hence, 
9u{x, y) — cos ipj can be zero. Since cos 9v{x, y) =r cos(</),y — <j)) is bounded by both 
r < cos{'tjji — Em) and cos(^y — £„) as |(^| < it follows that 

Oi.{x, y) > max{V';, (pi.} - Em ■= zi^^ - > 0. 
For each u, we choose an index such that 

\0„{x,y)-ipj^ \ < 
We need an estimate on jV. Since 

2 sin^ '^jV ^ _-|^ _ QQg^^^.^ _ g-^^ < 1 — cos(^^(a;, j/)) = 1 — r cos{(pv — (p) 

= 1 - r + r(l - cos{(pi, - (p)) <2 sin^ y + 2 sin^ 
we obtain {(pj^ — Em)"^ < ^(''Pf + ('/'i' ~ from which it follows that 
(5.12) j.-l< IVP + J2^ < ^{1 + 2u). 
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Using the fact that for j ^ j^, 
we obtain the estimate 

^ m/2 rn 

= 2m + 1 ^ ^ ™' 



m/2 . 2 , 



sin((2m+l)(6l^(x,y)-7/;j)) 



4 sin 0i/(a;, j/) sin sin 



Asiii9„{x,y) sin 



2 

1 \ ^ ■^-^ sm tpj 



' + 1 § 4 sin y) sin ''-^"f"^^' sin M£:|)±^ 

2 ; 2 1 / 

~ 8 - ern 2(2m + 1) - £m ^ {-ipj^ - V'il - 

Using the inequahty (|5.12l) . the first sum in the last expression is bounded by 



Q m/2 ^ „ 



-y 



< — m, 



< 



16 max{/, 2iy}-l/2 - 8 
while the second sum is bounded by, again using H5.12|l . 
2 '"/^ -, /i^-i . m 

^ y ^ y 1 + y _ 

2(2m + 1) ;^ Zi., ~ e„, \^ ^ - j - 1/2 ^^^^ j - - 1/2 ^ 

m/2 ^ I m-j„ . ^ . 

1 g max{/,2.}-l/2 [^^ ''^^'^ + + g T^2 ^ 

^ ? II SI o X 1/0 (^'"^ l^S^-^'" + 1) + ("1 - + (> + 1/2) log(m - > + 1)) 

2 maxjt , 2i/| — 1/2 

m/2 m/2 

< -(tT + l)10g(m + 1) > r -+?71> r 

- 2^ ^ max{?,2iy} - 1/2 max{/, 2;^} - 1/2 

< (tt^ + TT + l))mlog(m + 1). 

This completes the estimate for L^i\- Similarly, but more easily, we have by H5.12|l . 

m/2 . 2 , , ™/2 , 

" h 4sin2 Mf^l^ + 4(2m + 1)^ A. ^2^^ ^.^2 

TO TO x ^ 1 TO TO x ^ 1 ™ / 1 

This competes the estimate of 

Putting all cases together completes the proof of the lemma. It is easy to see 
that the largest constants in front of to log(TO + 1) or to in the three cases come 
from Case 3. □ 
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